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1 Abstract 

The authors used the embedded boundary method (EBM) to solve the two- 
phase incompressible flow with piecewise constant density. The front tracking 
method is used to track the interface. The fractional step methods are used 
to solve the incompressible Navier-Stokes equations while the EBM is used in 
the projection step to solve an elliptic interface problem for the pressure with a 
jump equal to the surface tension force across the interface. Several examples 
are used to verify the accuracy of the method. 

2 Introduction 

We consider here the two-phase incompressible Navier-Stokes equations (Zhijun 
Tan et al. gU) 

p{ut + {u-V)u) = -Vp + V ■ n{Vu + {Vuf)+F{x,t)+g{x,t), (1) 
V-w = 0, 

where p is the density (constant in each phase), u is the velocity field, p is the 
pressure, /i is the viscosity, g{x,t) is the external body force and 

F{x,t)= f f{s,t)S{x-X{s,t))ds 

is the singular interface force, concentrated on the interface with parametrization 
X{s,t). 

There are many different methods for solving the two-phase incompressible 
flow that also track the phase interface. Some of the most popular methods are 
the immersed boundary method [27], the front tracking method [37j .35J, the 
level set method [331 l30l[26] . and the volume-of- fluid method [28l[T5]. To deal 
with the singular source term F{x,t) such as the surface tension in equation 



ID) , these methods simply use a discrete delta function to transform the source 
term F{x,t) defined only on the interface to the grid cells near the interface, 
and then solve the Navier-Stokes equations as a one-phase problem. 

Recently, sharp interface methods to solve the two-phase flows directly us- 
ing the jump conditions due to the singular force term have received attention. 
The ghost-fluid method |12[ [23J [TB] has been used to solve the elliptic bound- 
ary value /elliptic interface problem, and then used to solve the two-phase in- 
compressible flow to more accurately compute the solution satisfying the jump 
conditions for the Navier-Stokes equations across the material interface. How- 
ever, since the ghost-fluid method uses approximate jump conditions, it can only 
achieve at most first order accuracy for the elliptic interface problem. Another 
popular method is the immersed interface method [121 [201 [HI [M] which has 
been applied to solve elliptic interface problems, parabolic interface problems, 
and two-phase incompressible flow with discontinuous viscosity. However, it 
seems that this method has not been used to solve variable density two-phase 
incompressible flow. And in addition to the basic jump conditions, the immersed 
interface method requires higher order derivative jump conditions which are of- 
ten difficult to be derived. 

In this paper, we use the embedded boundary method (EBM) [13] to solve 
the two-phase incompressible Navier-Stokes equations with piecewise constant 
density. This method was originally proposed [13] to solve the elliptic boundary 
value problem on irregular domains. It has since been used to solve the heat 
equation [^ [5Sj and the incompressible flow on a time-dependent domain [2S] 
with second order accuracy. Recently it has been extended to solve the elliptic 
interface problem with second order accuracy in two and three dimension |39| 
[^ where the jump conditions for the potential and its flux are used in the 
discretization of the elliptic equations using the finite volume method. 

The front tracking code FronTier developed at Stony Brook University 
has been successfully used for solving compressible flow |11[ [9l [5] . It uses the 
ghost-fluid method [121 \n[ [T5] to obtain sharp interface solution. In this pa- 
per, we solve the two-phase incompressible flow equations using the embedded 
boundary method, with the front tracking method used to track the material 
interface. The jump conditions due to the surface tension force are accurately 
solved in the projection step using the EBM while the jump conditions due to 
the discontinuous viscosity are disregarded for simplicity. The main contribu- 
tion of this paper is the development of a sharp interface EBM for the two phase 
incompressible flow. We also demonstrate the extension of the method to the 
cylindrical coordinate. An extension of the method has been used to solve large 
density ratio incompressible multiphase magnetohydrodynamic flows [T]. 

The paper is organized as follows. In section [31 we review briefly the pro- 
jection method used to solve the incompressible flow without considering the 
interface. In section [H we describe the embedded boundary method used to 
solve the incompressible flow with an internal interface. In section [SI we show 
some examples to verify our methods. Finally, we give the conclusion. In the 
appendices, we discuss the consequence of not considering the jump condition 
due to discontinuous viscosity, the code structure, and parallelization. 



3 The Projection Method 

In this section, we give the basic description of the projection methods used 
to solve the incompressible flow without considering the interface. Zhou et al. 
[10] contains the verification of the basic code for solving one phase flow and 
an implementation of the two phase flow using the immersed boundary method 
with the front tracking method. 

To solve the Navier-Stokes equation, one of the most popular methods is 
the projection method. There are many variations of the projection method 
[6]. In this paper, we use two different but similar projection methods (denoted 
as PMl and PM2 later on) with the embedded boundary method to solve the 
two-phase incompressible flow (PMl is first order accurate and PM2 is second 
order accurate when there is no interface). The projection methods for solving 
the Navier-Stokes equations ([Ij consists of two steps [6,. One first solves the 
convection and diffusion term to compute an intermediate velocity at the new 
time step. One then solves the Poisson equation to obtain the pressure, and 
use the new pressure to calculate the new time step velocity which satisfies the 
divergence free condition. One salient feature of the two projection methods 
used in this paper is that the pressure instead of the pressure increment is 
calculated in the projection step. This feature makes it easier to use the EBM. 
The reason is that for two phase incompressible fiow with surface tension, it is 
the pressure (not the pressure increment) which has a jump equal to the surface 
tension across the interface. 

3.1 Projection Methods 

We describe three different variations of the projection methods. However, the 
last projection method (PM3) is presented here only for comparison purpose. 
For simplicity, we use L to denote the discretization of the heat operator in the 
Navier-Stokes equations. 

3.1.1 Projection Method PMl 

This projection method is similar to Chorin's original projection method [71 [B]. 

u* — L{u", —u ■ Vu) 

7 = — Vp"+5 (2) 

At p ^ ^ 

where n is the time step index. At is the time step size, u* is the intermediate 
velocity, and u"+^ is the velocity at the new time step satisfying the divergence 
free property. Note that this method is only first order accurate. A notable 
characteristics of this method is that the projection step solves for the pressure 
p"+2 instead of the pressure increment. 
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3.1.2 Projection Method PM2 

A second order accurate method [53] is the foUowing 

(3) 
(4) 

(5) 
where the projection step solves for the pressure p"+2 (the same as PMl). 

3.1.3 Projection Method PM3 

Another popular second order accurate method is given by Bell, et al |310: 



= -iv0"+3 (6) 



At p 

Note that the projection step solves for the pressure increment 0"+2. For this 
reason, it is difficult to be used with a jump condition for the pressure itself 
(such as the surface tension between two phases). Thus, this method is not 
used in this paper. 

3.2 Time Discretization for the Calculation of the Inter- 
mediate Velocity 

There are many different methods to discretize the heat operator. In the fol- 
lowing, we briefly describe the Crank-Nicolson method and the Additive Runge- 
Kutta Method that we have used in our code. 

3.2.1 Crank-Nicolson Method 

Using the Crank-Nicolson method for the heat operator and using PMl as ex- 
ample, we have 

/ II* — 7;"\ 1 1 

where the convection term is solved using an explicit Godunov-type scheme as 

in 141: 



Table 1: Additive Runge-Kutta coefficients tableau 
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extrapolate the velocity to the cell face at time n+ ^ (the following shows 
the formula for only one face): 
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(8) 



where i, j denote the cell index, (i + 1/2, j) denotes one of the the cell face, 
n denote the time step index, and % can be replaced by the Navier-Stokes 
equations 

ut = -(u ■ V)m Vp + /i(V^u) + g(x, tj). 

P 

• solve the Riemann problem (Burgers' equation) to find the cell face veloc- 
ity at time n + 5. 

• use the cell face velocity to calculate the convection term (u ■ Vu)"+2. 



3.2.2 Additive Runge-Kutta Method 

The implicit Runge-Kutta method ([MJ) was used to solve the time dependent 
parabolic initial boundary value problem in (jlll[29j) and the parabolic interface 
problem in [35]. In this paper, we instead use the additive Runge-Kutta method 
[21| which seems to be easier for the discretization of the Navier-Stokes equation. 
Using the notation in [21j , we want to solve the following ordinary differential 
equation (ODE): 

V'{t) = l{t,y)+g{t,y) (9) 

where / is linear operator of y and 5 is a nonlinear operator. We use an implicit 
scheme for / and an explicit scheme for g. When used with PM2, / = Au 
and g — —u ■ Wu — Wp"^^ , Note that to solve the convection term u ■ Vu for 
the additive Runge-Kutta method, there is no need to do time extrapolation in 
equation ([8]). 

To solve the ODE ([9]), the s-stage Runge-Kutta method has the following 
form 

s s 

vt^ = 2/i""'^ + /^E «'J-/(^"-i + ^J-'^' yf^^ + hJ2b,,g{t^^i + c,h, yf) 



j=i 



j=i 



where i = 1,2, ...,s, n = 1,2,3,..., and Ci = X]5=i '^jj — X]^=i^y- The coeffi- 
cients are generally written as the tableau in Table [T] 

In this paper, we use the L-stable two stage additive Runge-Kutta scheme 
shown in Table [31 



Table 2: L-stable two stage additive scheme 
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3.3 Calculation of the New Velocity 

Now we consider the projection step of PMl and PM2 for calculating the pres- 
sure and the new divergence free velocity. For a regular grid cell containing no 
interface, the new velocity is calculated using: 



-Vp"+i/2 (10) 



At 



where p"+i/2 denotes the new pressure and is calculated using 
and the boundary condition 



(11) 



n ■ Vp"+i/2 ^ Q^ (12) 

We have used V • u"+^ = and equation (ITUl) to obtain the Poisson equation 
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4 The EBM for Two-Phase Incompressible Flow 

In this section, we first review the jump conditions of the incompressible flow 
with an internal interface. Next we review our embedded boundary method for 
solving the elliptic interface problem in subsection 14.21 Then we show how to 
use the EBM in a cylindrical coordinate system. The method to calculate the 
pressure gradient used in PM2 is described in subsection 14.41 At last we apply 
the EBM to solve the elliptic interface problem of the projection step for the 
two-phase incompressible flow in subsection 14.51 

4.1 The Jump Conditions across an Internal Interface 

For two phase flow with discontinuous coefficients, the solutions satisfy the 
following jump conditions across the internal interface: 







(13) 



du 



(14) 



[pJ.n^O (16) 

[-|^]-0 (17) 

where n is the normal vector and r is the tangential vector. Note that we have 
4 jump conditions for the two velocity components in 2D (6 in 3D) and 2 jump 
conditions for the pressure. 

• Jump condition (ITS)) is due to the viscous flow (Zhi, et al [34] ) . 

• Jump conditions ([T^ [TS)) are the result of force balancing in the normal 
and tangential directions (Ito and Li |13|1. The derivation of the jump 
conditions (|14ll5p is given by Ito and Li [T3] for the Stokes equations. It 
is trivial to extend it for Navier-Stokes equation. 

• Jump condition (ITS)) is due to [V • u] = (Lai and Li \n\). 

• Jump condition ((TTj) is needed to obtain a continuous velocity in the nor- 
mal direction across the interface when solving the elliptic interface prob- 
lem for the pressure. 

In this paper, we assume that the singular force term on the interface has only 
a normal component, which means that / ■ t = 0. The surface tension force 
satisfies this assumption. For simplicity, we also ignore the jump conditions due 
to the discontinuous viscosity. Appendix [Xj shows the consequence of such a 
simplification. Therefore the velocity and its derivatives are continuous across 
the interface, while the pressure has a jump equal to the surface tension across 
the interface. 

4.2 The Embedded Boundary Method for the EUiptic In- 
terface Problem 

We review briefly the embedded boundary method for solving the elliptic in- 
terface problem. For more detail, refer to ISHl ISl [IH HI] • We assume that the 
interface can only cross any cell edge at most once, as implicitly assumed in 
the Marching Cubes algorithm |23j. This assumption is necessary to limit the 
number of different cases that could arise. The geometric information needed of 
the partial cells needed by the embedded boundary method is calculated using 
the divergence theorem for each cell case by case |38) . 

The elliptic interface problem is a special elliptic problem with an internal 
interface: 

V.^^/ (18) 



Table 3: Number of material components and cell unknowns for different cell 
types 
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Components 
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External 


N/A 


N/A 


N/A 


Internal 


1 


1 





Boundary 


1 


1 





Partial 


2 


2 


2 



where p is a piecewise continuous function with jump across the internal interface 
and / is a given function which is continuous inside each part of the domain. To 
close the problem, boundary conditions are needed for the exterior and interior 
boundary. Either Dirichlet or Neumann boundary can be given on the exterior 
boundary. For the interior interface, we have the following two jump conditions: 



[p] = Ji(x), 
p on 



(19) 
(20) 



where Ji and J2 are given functions of the spatial variables [T9]. Note that our 
method could be easily extended to more general cases where Ji and J2 are 
functions of the unknown variables p defined on the interface. 

The EBM can be used to solve the elliptic problem with an irregular domain 
boundary and an internal interface. It uses a Cartesian mesh. The mesh cells 
are classified into four types: 

• An external cell is outside of the computational domain and thus is not 
used in the computation. 

• An internal cell is a cell wholly located inside the computational domain, 
possibly with one of its cell face being the domain boundary. 

• A boundary cell is a cell intersected by the irregular exterior domain 
boundary with part of the cell out of the domain and part of the cell 
inside the domain. 

• A partial cell is a cell intersected by the internal interface. It is separated 
into two or more parts by the internal interface. Those different parts are 
also called partial cells in the following. 

When the EBM is used to solve the elliptic interface problem, one or more 
unknowns are defined at the cell center, as shown in Table [3] For an interior cell 
or a boundary cell, only one unknown is needed at the cell center and a standard 
finite volume method can be used to setup one algebraic equation using the 
elliptic equation as in j39[ [T4l [24l [29] . For a partial cell, four unknowns in total 
are needed to make the discretization of the elliptic equation consistent with 




Figure 1: Placement of unknowns in a partial cell containing the cell internal 
interface, pa and pb are cell center unknowns for component a and 6 respectively. 
Similarly, Pintfc.a and Pintf,b are cell interface unknowns for component a and 
h respectively. 

the two interface jump conditions ([T^ . (I^Ul) . Figure [T] depicts the placement of 
unknowns in a partial cell with internal interface. The cell contains two partial 
cells (for material components a and b) which are separated by the interface. 
For each partial cell, there is one unknown defined at the cell center {pa for 
component a, pf, for component b) for the discretization of the elliptic equation. 
In order to satisfy the two jump conditions (jl9[|20p . two more unknowns Pintfc,a 
and Pintfc,b are defined at the center of the cell internal interface (portion of the 
interface contained within the cell) for the two components a and b. Thus, four 
unknowns are defined in total for one partial cell. And four algebraic equations 
can be constructed using the elliptic equation for the two partial cells and the 
two jump conditions across the cell internal interface 

4.2.1 Discretization of the Jump Conditions 

We first describe the method for the discretization of the two jump conditions 
across the interface for the cell (i, j). A schematic of the corresponding stencil 
and states used in the interpolation method is shown in Figure[5] Two unknowns 
Pintfc,a and Pintfc,b are defined at the center of the cell interface for components 
a and b respectively. 

We assume the direction of the normal to the interface as pointing from a 
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Figure 2: Stencil for the discretization for a partial cell 



to h. For the first jump condition (|19p . the discretization is simply 

Pintfc,h Pintfc.a — <-'l 



(21) 



To discretize the second jump condition (|20p . we need to calculate the normal 
direction derivatives of the unknowns. There are many different approach avail- 
able as in [39l [Ml [24] . The main idea is to construct a polynomial using the 
unknowns with the same component and then take derivative along the normal 
direction to get the flux for that component. To construct the polynomial for 
component a, we pick cell interface unknown Pintfc.a from cell (i,j) and then 
cell center unknowns from other neighboring cells. In 21?, we need 3 unknowns 
to construct a linear polynomial (for first order accurate flux) and 6 unknowns 
to construct a quadratic polynomial (for second order accuracy flux). Least 
square fitting is also possible. 

Taking the normal derivative of the constructed polynomials, we obtain the 



normal derivatives at the cell interface center for the component a, -J^ 



the component 6, 



dp 
dn 



Using the second jump condition (|20l) . we have 



dp 
dn 



dp 
dn 



and 



(22) 



Method for calculating the flux across the cell internal interface We 

use linear polynomial construction in 2D as an example. A linear polynomial 
in 2D can be written as 

p{x, y) ^ l3o + Pix + Ii2y 

where /3o, Pi and /32 are undetermined coefficients. Thus we need to find 3 
cell unknowns (1 unknown should be a cell interface unknown) for constructing 



11 



the polynomial. Denoting the coordinate locations of the unknowns as {xi,yi), 
(3^2,2/2)7 and (x3, j/3), the unknowns as pi, p2, and ps, we have 




Denoting the above matrix as X, the coefficient vector as /3, and the right 
hand side as P, we have 

or 

p = x-^p. 

Now the constructed polynomial can be written as 





The derivatives of the polynomial can be written as 

Px = Pi 
Py == ^2 

— = p^Tii + pyn2 = ( ni n2 ) ( ^1 ) = n^/3 

where n = ( ni, n2 )^ is the unit normal vector, and h = [ 0, ni, n2 ) ■ 
Thus we have 

1^ = n^X-^P 
on 

= {{x-Yufp, 

= {{X^r^nfP 

Therefore, we can represent g^ as a linear combination of the unknowns P. Sim- 
ilarly, we can calculate the directional derivative using second degree polynomial 
construction. 

Now suppose that instead we want to use least square fitting, we have 

(3 = {{X^X)-'X^)P 

The directional derivative can be calculated as 

^ = h^(3 = n'^{{X^X)-'X^)P 
on 

= (xix^xy'hfp, 

thus 3^ can also be written as a linear combination of the unknowns P. 

on 
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4.2.2 Equations for the Cell Center Unknowns 

For the two cell center unknowns defined at the partial cell {i,j), we can use 
the EBM technique to set up two algebraic equations by integrating the elliptic 
equation over the corresponding partial cells. See Figure [5] for the stencil to set 
up the equation for the unknown of the component a using the partial cell cdef. 
Integrating the equation (J18p over the partial cell cdef and using the diver- 
gence theorem, we obtain the following expressions: 

Vp I ^P ff 

V • — dxdy ~ (p • nds ~ fdxdy, 

cdef P Ja(cdef) P J J cdef 

or 

• nds + / • nds + / • nds + / • nds ~ fdxdy, 

d P Jde P Jef P Jfc P J J cdef 

which is 

led ■ fluXcd + Ide ■ fluXde + lef ' fluXef + Ife ' fluXfc = / / fdxdy (23) 

J J cdef 

where Imn is the length between m and n. Therefore we only need to calculate 
the flux across the cell edges or cell interface. For fluXcd, a second order deriva- 
tive is calculated by using a linear interpolation of ^''^"^~^''^ and ^'^^'^~^~^'^^'^ 
to the center of cd (see |391I14|). For fluxde, we simply use P'+i^ P''J to calculate 
the derivative. For fluxef, a linear interpolation of ^''^^^'^''^ and p^+^-^+J^^p^+^-j 

to the center of ef is used. And fluxfc is calculated by ^|a used in (P^ . Note 
that we need to multiply the derivatives calculated with - at the corresponding 
edge or cell interface center to obtain the flux. In the same way, we calculate 
fluxes for the other partial cells. More details can be found in [39lll4| . 

4.2.3 Geometric Property Calculation 

We use 3D as example. Since the EBM is a finite volume method, we need 
to calculate the partial cell volumes accurately. To solve the elliptic interface 
problem, we also need to calculate the cell interface area, normal and center. 

Due to the assumption that the interface crosses the cell edge at most once, 
there are finite number of cases for the partial cell configurations with two 
components. The cell volume, cell interface area, interface normal, and interface 
center are calculated case by case by using the divergence theorem: 



V • Fdv = / It ■ Fds. (24) 

n Jan 

where F = {Fx,Fy,Fz). Thus, we can use a cell boundary integration instead 
of a cell volume integration to calculate the geometric information needed. 
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In the Cartesian coordinate, the divergence operator is defined as 

dx dy dz 

To calculate the volume of the domain 17, we can let -F = X = (a;, y, zY' 
where X is the Cartesian coordinate. Note that F is not unique. Then we have 

V • Fdv = In- Fds 
n Jon 



Jn Jon 



m Jdn 

Therefore, we can calculate the volume using the surface integration 

f r„o 11 ■ Xds , . 

Volume ^ dv^ JdO__^ _ (25) 

Jn 3 

For the cell interface area, normal and center, we do surface integrations: 

ds 

intfc 



n ds 

intfc 

I ~ids 

J intfc 



where intfc refers to the cell interface. 
For more detail, refer to |38j . 

4.3 The Embedded Boundary Method in Cyhndrical Co- 
ordinates 

In this subsection, we describe the embedded boundary method used to solve 
an elliptic interface problem in a cylindrical coordinate. 

To discretize the elliptic equation (|18|) for a mesh cell, we use the divergence 
theorem to get 

® • nds = / fdv 

or 

---ds = / fdv 
p on J 

This is true for all coordinate systems. However, we need to change the 
formula for the calculation of the geometric property for the mesh cell, such 
as the cell face area, cell volume, cell interface area, normal and the center. 
The formula for the flux calculations across the cell face and the cell internal 
interface also need to be modified. 
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4.3.1 Finite Difference Scheme for the Flux Calculation 

To calculate the flux across the cell faces, we use 

dp _ 
on 

where the gradient operator is defined as 

or r ou oz 

and n is the normal of the cell faces. 

For the flux across the partial cell interface, one way is to use a modified 
version of the method presented in subsection 14.21 for Cartesian coordinates. 
Another first order method is the following. We have 

dp = Vp • nds. 

Therefore 

Ap w Vp • nAs. 

Thus we have 

Vp • n sa — — 
As 

4.3.2 Geometric Property Calculation 

The surface elements in cylindrical coordinate for the cell face of the internal 
cell are 

ds = rdOdz, 

ds = drdz, 

ds — rdrdO 

in the respective coordinate planes. They can be used to calculate the cell face 
area for the 6 cell faces. The volume element is 

dv — rdrdOdz. 

To calculate the partial cell volume, we also use the divergence theorem. In 
a cylindrical coordinate, the divergence operator is defined as 

r dr r 86 dz 

To calculate the volume of the domain fl, we can let F — {r,r9,z)'^, and 
then we have 



/ V • 'fdv = n- 'fds 

Jn Jdfi 

Adv = / '^ ■ ~^^^ 
Jdn 
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Note that the choice for F is not unique. Therefore, we can calculate the volume 
using the surface integration 

/■ Lo n ■ Fds , , 

Volume ^ dv^ JX]_ _ (26) 

For more detail, refer to |38| . 

4.4 Calculation of the Pressure Gradient 

To calculate the pressure gradient used in ^ and Q, we simply use least 
square fitting. To calculate the gradient for cell {i,j) with component c, we do 
the following: 

1. find the nearest cell {i,j) containing component c. 

2. collect all the states with the same component from neighboring cells with 



cell index {i,j) satisfying 



< 2 and 



J -J 



< 2. 



3. use least square fitting to fit a quadratic polynomial and take gradient to 
calculate the pressure gradient needed. See subsection 14.2.11 for detail. 

4.5 The EBM for the Projection Step 

We solve the elliptic equation (fTTj) in the projection step. The EBM for a 
general elliptic interface problem is given in previous subsections. The two 
jump conditions given by P3 [T7)) are used to connect the pressure solution 
across the interface. 

Thus, for the projection step, the equation is 

V . (iv(p"+i/2)) - ^V • w* (27) 

coupled with the two jump conditions (omitting the effect due to the viscosity 
jump across the interface) 

[p"+i/2] =f-n (28) 

and 

where f ■ n ~ an for the surface tension force between the two phases {a is the 
surface tension coefficient and k is the mean curvature). 
Our projection method consists of the following steps: 

1. Calculate cell face velocity uZ^^^ using intermediate cell center velocity u* 
by linear interpolation. 
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2. Calculate the elliptic interface equation for p: 

^ At^ ■ '"■*face 




where the right hand side of the elliptic equation is calculated using the 
divergence theorem. 

3. Calculate the projected cell face velocity u^'^^^ ^ for each full/partial cell 
face using: 

,„+i ,.* At dp 



u 



f""''-'' ^"^^ Pface.c ■ dn 



face^c 



where S/pface,c is calculated using the cell center pressure for component 
c. 

4. Calculate the averaged cell face velocity u'l^^^ using u^^^^ ^. For example, 
in Figure m 

"se — , 

l-ge 

where Ige-, Ifg, Ife are the length between corresponding point on the cell 
face, u^tl is the velocity through gf and u'l'^]^ is the velocity through /e. 

5. Calculate the cell centered velocity using u'i^^g by using the 2nd order 
TVD reconstruction algorithms [2]. A simple alternative is to use inter- 
polation of the cell face velocity to calculate the cell center velocity. 



5 Examples 

We compare theoretical and simulated bubble oscillation frequencies in 2D and 
2iD to verify our embedded boundary method to solve two-phase incompressible 
flow. To apply EBM to cylindrical coordinates, we first verify the accuracy of 
the EBM by solving an elliptic interface problem with known solution. Then we 
apply our methods to solve a more complicated problem of engineering interest. 

5.1 Bubble Oscillation in Two and Three Dimension 

We consider the oscillation period of a droplet under zero gravity, for which 
there is an analytical solution and the surface tension is the dominant force. 
We have used both PMl and PM2 to obtain the simulation result. However, it 
is found that the difference between the results obtained using the two methods 
are negligible. 
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Table 4: Bubble Oscillation in 2D, domain [0,2] x [0,2], p 
{0.0005, 0.0005 X 0.01}, a = 0.5, Rq = 0.8, n = 2 and e = 0.05. 
period Tz = 2.56638189 



= {1,0.05}, V = 
with theoretical 



Mesh Size 


Period 


Error in Percentage 


20x20 


3.4602 


17.45 


40x40 


3.0476 


9.95 


80x80 


2.8962 


6.38 


160x160 


2.7643 


3.93 



5.1.1 Bubble Oscillation in Two Dimension 

For a 2D droplet under zero gravity, when the initial position of the droplet 
interface with small perturbation is given by 

R{e) = Ro + €cos{ne) 

where Rq the unperturbed radius, e is the amplitude of the perturbation and 
n is the order of the Legendre polynomial, the oscillation frequency is given in 
[inilS] as 

/ [n^ — n)cr 



y {pd + Po)Rt 

where pd and po are the densities for the droplet and outer fluid and a is the 
surface tension coefficient. 
The period is given by 

2tt 

We run the simulation with the domain as [0,2] x [0,2], the densities and dy- 
namic viscosities of the droplet and the outer fluid given hy p = {1,0.05}, 
z/ — {0.0005, 0.0005 X 0.01}. We impose a surface tension coefficient of cr = 0.5. 
For the initial interfacial position, we use Rq = 0.8, n = 2 and e = 0.05. Table|4] 
shows the convergence of the oscillation period under mesh refinement. Figure[3] 
shows the interface velocity of the droplet at time t — 0.5. From the picture we 
can see that the interface velocity changes smoothly along the interface. Figure 
|4] shows the pressure over the whole computational domain. It is apparent that 
the pressure has a jump across the interface due to the surface tension. Due 
to the small perturbation of the initial interface (e = 0.05) and large radius of 
the droplet {Rq — 0.8), the variation of the surface tension along the interface 
is very small. Therefore, the variation of the pressure jump along the interface 
is small and not apparent in the figure. We observe the accurate resolution of 
the pressure discontinuity up to the interface without over shoots, the Gibbs 
phenomenon. 
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Figure 3: Bubble Oscillation in 2D, with interface velocity drawn as vector 
starting from the interface at time t — 0.5 




Figure 4: Bubble Oscillation in 2D, Pressure at time t = 0.5 
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Table 5: Bubble Oscillation in 3D, theoretical period is T2 — 2.222 



Mesh Size 


Period 


Error in Percentage 


20x20x20 


2.433 


9.5% 


40x40x40 


2.300 


3.5% 


80x80x80 


2.184 


1.71% 



Figure 5: Bubble Oscillation in 3-D, Velocity Field on a Slice through the Center 
of the Droplet at i = 1.8 




5.1.2 Bubble Oscillation in Three Dimension 

For a 31? droplet under zero-gravity, when the initial position of the droplet 
interface with small perturbation [32] is given by 

r{e,t)=Ro + ePn{cos{e)), 

where Rq is the radius of the droplet, P„ is the Legendre polynomial of order 
n, e is the amplitude of the perturbation, then the frequency of the droplet 
oscillation is given by 



UJ„ 



' 1 n{n - l){n + l){n + 2) 
We Rl{n + l + nX) 



— — is the Weber number. 

We perform the 3D simulation with the domain as [0, 3] x [0, 3] x [0, 3]. We 

use Re = ^^^ = 2000, We = l^^^^ = 1 and pa = 1, po = 0.001, p.d = 0.0005, 
Ho = 0.000005. Table [5] shows the convergence of the oscillation period under 
mesh refinement. Figure [5] shows the velocity field on a slice through the center 
of the droplet. 



where X — —, and We 
pi 
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Table 6: Mesh Convergence Study for the Elhptic Interface Problem in Cylin- 
drical Coordinates 



Mesh Size 


Lao Error 


L2 Error 


Li Error 


10x10x10 


0.03376388 


0.00184484 


0.00056186 


20x20x20 


0.01980214 


0.00107820 


0.00034137 


40x40x40 


0.01104179 


0.00059434 


0.00019211 



5.2 EBM for the Elliptic Interface Problem in Cylindrical 
Coordinates 

We use the method of manufactured solutions to verify our EBM implementation 
for the elliptic interface problem in the cylindrical coordinates. The computa- 
tional domain is r e [1,1.628], 9 € [0,0.628] and z S [0.628]. The interface 
position is a sphere relative to the cylindrical coordinates, given as 



v/(r - 1.314)2 + (6» - 0.314)2 + (^^ _ 0.314)2 = o.2. 

1.03 for the density 



/ 



(30) 



We use pd — 0.811 for the density inside the sphere and po 
out side. We solve the equation 

P 
where / is a given function. We use 

p(r, tl, z) = e 5 

as the exact solution and substitute into the elliptic equation ([50)1 to calculate 
the right hand side /. Table [5] shows the mesh convergence for the this problem. 
From the table, we can see that the method is first order accurate, due to our 
use of a first order method to calculate the flux across the interface. 

5.3 High Speed Two Phase Couette Flow 

Here we use the EBM to simulate high speed two-phase Couette mixing in a 
31? angular sector. A more detailed description of this simulation using the 
immersed boundary method can be found in |40j . In the current simulation, 
we use paqu — 1.03g/cm^ and porg = 0.811g/cm^ for the fluid densities for 
the aqueous and organic phase respectively, paqu = 0.0l02g/cm ■ s and porg = 
0.016g/cm ■ s for the viscosities, and a = lOdyn/cm for the surface tension. The 
computational domain is an angular sector of the cylinder with r e [2.538, 3.166], 
e [0,0.314], z £ [0,0.628]. Figure [H shows the interface position at t = 28/zs. 

The explanation of the result will be given in another paper. This example is 
given here to show the capability of the implemented EBM to deal with problem 
with complex interface and of engineering interest. Under the assumption that 
there is at most one interface crossing on each cell edge, there are 2^ — 256 
possible cases for a mesh cell. The numerical simulation results for this example 
show that all 256 cases appeared at the same time. 
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Figure 6: Interface for Two Phase Couette Flow 




6 Conclusions 

In this paper, we extended the embedded boundary method to solve the two- 
phase incompressible flow. We verify our method by solving the droplet oscilla- 
tion problems in 2D and 3-D. We also show that the EBM can be easily extended 
to solve the elliptic interface problem in other coordinate systems such as the 
cylindrical coordinates. Finally, we simulated the interface contact problem in 
a rotating cylinder to show the robustness of the method. Currently, we did 
not consider the jump condition due to the discontinuous viscosity. This will be 
addressed in the future. 
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A Approximate Jump Conditions for Decoupling 
Velocity 

This part is modified from [TB]. Instead of using the jump condition (fTO|) . we 
can use any equivalent jump conditions. Due to the incompressibility condition 
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V • w = 0, we have [/iV • u] = 0, which can also be written as (Zhi, et al |34) ) 
r du. , du . r du . ^ .„^, 

Equations (I15l3ip can be written in matrix notation in 3D as the foUowing: 

du , 




'^^ ^ du 





Multiplying both sides of this equation from the left by 

T 



we have 







dn 








Since the tangential derivative of the velocity is continuous ([17], [K]), we have 




T / \ / \ T 

Tl 



>"£> -"'>;■ (^l^M") i"]^, '^'' 








0(34) 



Thus if the surface force / has no or very small tangential components (as 
in surface tension force) and [fj] is very small, it is safe to approximate the jump 
conditions of the velocity using 

[^^^] = (35) 

At an interior point away from the interface, we have [^] = and / = 0, 
therefore [/ifi] = 0. 

B Code Structure and Parallelization 

In this section, we briefly introduce the data structure, coding for dealing with 
partial cells, and the parallelization method used for the EBM. For simplicity, 
we use 2D as example. 
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The EBM uses a structured Cartesian mesh. For a 2D rectangle domain, we 
use a two dimensional array to represent the Cartesian mesh. We need to store 
extra geometric information for partial cells (partial area, cell interface length, 
center, normal) and partial cell edge (length, center). We also need to store extra 
states for the partial cells and partial cell edges. Since the data size needed for 
partial cells with internal interface is much larger than the data size needed for 
a regular internal cell, we use an extra pointer to a data structure for the partial 
cells (called partial cell data structure later on). The extra data structure is 
allocated and deallocated dynamically. A cell type variable is maintained for 
each cell to distinguish between internal, partial, external and boundary cell 
(cell with external boundary crossing through). For a moving interface, the cell 
type could change dynamically. For the cell edge, the same data structure is 
used for both whole and partial cell edges since a partial cell edge needs only a 
few extra data storage. 

Whenever the interface is changed, the partial cells are identified and the 
partial cell data structures are set. The partial cell geometric information is 
calculated using cell corner component and cell edge crossing information using 
the Marching Cubes algorithm |231 [38] (with the assumption that there is at 
most one crossing at each cell edge). The geometric information for the partial 
cell edge is also set accordingly. 

The EBM is a finite volume method. A systematic way of enumerating the 
cell internal interface and the partial cell edges is used to calculate the flux of 
the differential equations. 

Due to the structured Cartesian mesh used, it is very easy to parallelize 
the code. There are three types of variables: cell center variables (cell center 
velocities, pressure), cell face variables (cell face velocities), and cell corner 
variables. The parallelization consists of three steps: 1) pack; 2) send/receive; 
3) unpack. The unpack step is just the reverse step of the pack step. We first 
allocate a big array. Then for each cell to be sent, we first pack the cell type, 
and then pack all other variables associated with that cell into the array. Then 
we send/receive the single array. The unpack process is similar with the pack 
process. For each cell in the buffer zone, we first get the cell type from the 
received array, and then unpack the other variables. 
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